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Executive  Summary 

We  report  here  on  the  work  performed  during  the  third  year  (April  2015  -  March 
2016)  of  the  contract  FA  8655-13-1-3061  on  Multiscale  materials  science:  a  math¬ 
ematical  approach  to  the  role  of  defects  and  uncertainty. 

The  bottom  line  of  our  work  is  to  develop  affordable  numerical  methods  in 
the  context  of  heterogeneous,  possibly  random,  materials. 

Many  partial  differential  equations  of  materials  science  indeed  involve  highly 
oscillatory  coefficients  and  thus  small  length-scales.  When  the  microstructure  of 
the  materials  is  periodic,  or  random  and  statistically  homogeneous,  homogeniza¬ 
tion  theory  can  be  used,  and  allows  to  appropriately  define  averaged  equations 
from  the  original  oscillatory  equations.  When  no  such  structural  assumption 
(e.g.  periodicity)  on  the  materials  microstructure  can  be  made,  homogenization 
theory  still  holds,  but  does  not  provide  any  explicit  formulae  amenable  (even 
possibly  after  some  approximation)  to  numerical  computation.  One  possibility 
is  then  to  directly  address  the  original  problem  (rather  than  passing  to  the  limit 
of  infinite  scale  separation),  and  to  use  dedicated  numerical  approaches  for  such 
multiscale  problems,  such  as  the  Multiscale  Finite  Element  Method  (MsFEM). 

In  this  report,  we  first  consider  a  multiscale  advection-diffusion  problem  on 
a  perforated  domain,  in  the  convection-dominated  regime  (see  Section  1).  There 
are  two  small  scales  in  the  problem.  The  first  small  scale  is  the  size  of  the 
perforations,  which  is  of  the  same  order  of  magnitude  as  the  distance  between  two 
neighboring  perforations.  The  second  small  scale  is  related  to  the  ratio  between 
the  diffusion  and  the  convection  coefficients.  In  addition  to  the  difficulty  owing  to 
the  presence  of  different  scales,  the  strong  convection  is  also  a  source  of  potential 
instabilities.  Taken  separately,  each  phenomenon  can  be  addressed  by  classical 
approaches:  MsFEM  type  approaches  and  stabilized  type  techniques  (e.g.  the 
Streamline  Upwind  Petrov- Galerkin  (SUPG)  method),  respectively.  We  show 
here  how  to  adapt  these  two  methods  in  a  unified  single  approach  to  efficiently 
solve  multiscale  advection-diffusion  problems  set  on  perforated  domains,  in  the 
convection-dominated  regime. 

The  problem  we  consider  here  is  similar  in  spirit  to  the  multiscale  advection- 
diffusion  problem  considered  in  [EOARD3,  EOARD6],  except  that  there,  the 
multiscale  nature  of  the  problem  was  stemming  from  the  fact  that  the  diffusion 
coefficient  in  the  problem  was  highly  oscillatory.  Here,  the  multiscale  nature  of 
the  problem  originates  from  the  geometry  of  the  domain  on  which  the  problem 
is  posed. 
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Recall  that  the  MsFEM  method  consists  of: 

•  an  “offline”  stage,  where  highly  oscillatory  basis  functions  are  numerically 
computed  as  solutions  of  local  problems  (that  mimick  in  a  suitable  way 
the  reference  problem  on  a  subdomain).  These  basis  functions  are  thus 
expected  to  be  well-adapted  to  the  problem  at  hand. 

•  an  “online”  stage,  where  a  Galerkin  approximation  of  the  reference  problem 
is  solved.  The  approximation  is  performed  in  the  finite  dimensional  space 
spanned  by  the  basis  functions  computed  in  the  prior  offline  stage. 

In  this  work,  we  compare  different  numerical  approaches: 

•  A  first  approach  (denoted  Adv-MsFEM  below)  is  to  take  into  account  both 
the  diffusion  term  and  the  convection  term  in  the  definition  of  the  basis 
functions  (see  (7)-(8)  below). 

•  A  second  approach  (denoted  MsFEM  below)  is  to  only  take  into  account 
the  diffusion  term  in  the  definition  of  the  basis  functions  (see  (10)-(11) 
below). 

In  the  online  stage,  both  approaches  may  be  stabilized,  or  not,  using  e.g.  the 
SUPG  technique  (see  (13)  and  (19)  below). 

Preliminary  results  have  been  obtained.  Our  partial  conclusions  are  as  fol¬ 
lows.  For  identical  discretization  parameters,  the  Adv-MsFEM  and  its  stabilized 
variant  (the  Stab- Adv-MsFEM  approach)  provide  results  with  similar  accuracy, 
which  are  better  than  those  obtained  using  the  MsFEM  or  the  Stab-MsFEM  ap¬ 
proaches.  The  costs  of  the  non-stabilized  variants  are  similar,  and  smaller  than 
the  costs  of  the  stabilized  variants,  due  to  the  need  for  the  latter  to  assemble  the 
stabilization  terms.  For  the  problem  of  interest  here,  the  method  of  choice  seems 
to  be  the  Adv-MsFEM  approach,  although  more  comprehensive  numerical  tests 
are  needed  to  confirm  this. 

This  ongoing  work  is  being  reported  on  in  [EOARD4].  It  is  performed  by 
Claude  Le  Bris  (PI),  Frederic  Legoll  (Co-PI)  and  Frangois  Macliot  (third  year 
PhD  student,  see  [EOARD7]),  thanks  to  the  specific  support  of  EOARD  to  the 
research  activity  of  our  group.  It  is  a  follow-up  of  the  work  [EOARD3],  described 
in  [EOARD6],  and  of  the  work  [EOARD2],  described  in  [EOARD 5]. 

We  next  describe  in  Section  2  some  works  related  to  the  construction  of 
coarse  approximations  for  problems  with  highly  oscillatory  coefficients.  Our  aim 
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is  to  define  and  construct  the  best  non-oscillating  coefficient  A  (think  e.g.  of  a 
constant  coefficient)  that  is  consistent,  in  a  sense  to  be  made  precise,  with  the 
behavior  of  a  heterogeneous  material  modelled  by  a  highly  oscillatory  coefficient 
A£(x).  Such  an  approach  can  be  considered  as  an  alternative  pathway  to  stan¬ 
dard  homogenization  techniques  when  these  latter  are  difficult  to  use  in  practice, 
in  particular  when  information  is  missing  on  the  coefficient  A£(x).  Of  course, 
modelling  the  material  with  A  (rather  than  A£(x ))  yields  much  more  affordable 
approaches,  since  there  is  no  fast  frequency  in  that  coefficient.  A  central  ques¬ 
tion,  of  course  intimately  related  to  homogenization  theory,  is  to  construct  in  an 
efficient  manner  this  best  constant  coefficient. 

We  focus  here  on  the  case  when  A£(x)  is  random,  which  is  a  discriminat¬ 
ing  case  where  the  standard  homogenization  approach  is  very  expensive.  We 
show  below  that,  for  an  essentially  identical  computational  cost  compared  to  the 
standard  homogenization  approach,  our  approach  allows  us  to  compute  a  more 
accurate  approximation  of  the  solution  u£  to  the  highly  oscillatory  equation  (both 
in  L 2  and  in  H 1  norms).  In  contrast,  it  is  less  efficient  for  the  approximation  of 
the  homogenized  matrix  A*. 

This  work  has  been  expanded  in  [EOARD1]  (it  is  also  briefly  described 
in  [EOARD6]).  It  has  been  performed  by  Claude  Le  Bris  (PI),  Frederic  Legoll 
(Co-PI)  and  Simon  Lemaire,  a  postdoc  student  fully  funded  by  this  contract  and 
who  joined  our  group  from  June  2014  to  September  2015. 
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1  Multiscale  advection-diffusion  problems  posed 
on  perforated  domains 

As  announced  in  the  previous  EOARD  report  [E0ARD6],  we  describe  here  on¬ 
going  works  concerned  with  multiscale  advection-diffusion  problems  posed  on 
perforated  domains,  in  the  convection-dominated  regime.  This  problems  brings 
together  two  difficulties  we  have  already  worked  on: 

•  problems  set  on  perforated  domains,  which  we  addressed  (in  a  purely  dif¬ 
fusive  context)  in  [EOARD2], 

•  advection-diffusion  problems  posed  on  standard  (i.e.  non-perforated)  do¬ 
mains,  which  we  addressed  during  the  PhD  thesis  [EOARD 7]  of  Frangois 
Madiot,  whose  research  activity  is  partly  funded  by  this  Grant. 

The  work  we  describe  here  is  a  follow-up  of  these  two  works,  with  a  strong 
practical  motivation. 

1.1  Introduction 

We  consider  a  bounded  domain  O  C  and  a  set  B£  of  perforations  within 
this  domain.  The  perforations  are  supposedly  small  and  in  extremely  large  a 
number  (so  that  they  are  close  to  each  other).  The  parameter  £  encodes  the 
typically  small  distance  between  the  perforations.  We  denote  by  Vt£  =  Vt\B£  the 
perforated  domain  (see  Figure  1). 


Figure  1:  The  domain  D  contains  (possibly  random)  perforations  Be.  The  perfo¬ 
rated  domain  is  Vt£  =  Vt\Be.  The  boundary  of  is  the  union  of  d  B£  D  1L  (the 
part  of  the  boundary  of  the  perforations  that  is  included  in  VLe)  and  of  <9DnDe. 
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The  problem  we  consider  reads 

—aA u£  +  I  .  =  /  in  f2e,  u£  =  0  on  dVt£.  (1) 

£ 

The  homogeneous  Dirichlet  boundary  condition  (i.e.  u£  =  0)  on  the  boundary 
of  the  perforations  is  an  important  feature  of  (1).  It  is  motivated  by  the  fact 
that  (1)  can  be  seen  as  a  step  toward  the  resolution  of  the  Stokes  problem  on 
perforated  domains.  In  that  case,  homogeneous  Dirichlet  boundary  conditions 
on  the  perforations  are  typical  for  many  applicative  contexts. 

The  parameter  a  in  (1)  is  a  positive  constant.  The  convection  field  b  is 
assumed  to  be  divergence- free.  In  the  sequel,  we  assume  that  the  ambient  di¬ 
mension  is  d  =  2  and  that  fl  =  (0,  l)2. 


Two  comments  are  in  order  regarding  the  definition  of  problem  (1). 

First,  note  that  we  have  scaled  the  convection  term  by  a  factor  1/e  in  (1). 
In  that  regime,  and  assuming  that  the  perforations  are  periodically  located  and 
that  b  is  periodic,  it  can  be  shown  that 


lim 

£ — >-0 


U£  —  £2  W 


f 


0, 


where  the  function  w  solves  the  so-called  corrector  problem 


(2) 


I  —a  Aw  +  b  ■  Viu  =  1  in  Y\B,  w  =  0  in  B, 
lye-)-  w(y)  is  Y'-periodic, 

where  Y  is  the  unit  cell  and  B  is  the  perforation  of  that  unit  cell.  Through  w,  the 
asymptotic  behavior  of  u£  hence  depends  on  the  convection  field  b.  In  contrast, 
if  we  consider  the  problem 


— aAv£  +  b{x/e )  •  'Vv£  =  f  in  v£  =  0  on  dfl£,  (3) 

where  the  convection  term  is  not  scaled  by  a  factor  1/e,  then  (2)  still  holds  but 
with  w  defined  by 

J—  a  Aw  =  1  in  Y  \  B,  w  =  0  in  B, 
lye-)-  w(y)  is  Y'-periodic. 

For  asymptotically  small  e,  the  solution  v£  to  (3)  can  hence  be  well  approximated 
by  a  quantity  independent  of  the  convection  field.  This  suggests  that,  for  e  small, 
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the  effect  of  the  convection  is  less  critical  (and  thus  less  challenging  to  capture) 
than  in  the  case  of  (1). 

Second,  we  recall  that,  in  [EOARD3]  (see  also  [EOARD6]),  we  have  consid¬ 
ered  the  multiscale  advection-diffusion  problem 

—  div  [Ae{x)'Wu£]  +  b{x)  ■  Vm£  =  /  in  hi,  u£  =  0  on  dQ,  (4) 

where  the  diffusion  coefficient  A£  varies  at  the  small  scale  £,  and  where  the  con¬ 
vection  vector  held  b(x)  is  given  and  large.  In  some  loose  sense,  the  consideration 
of  (4)  is  a  preparatory  step  toward  that  of  (1).  For  (4),  boundary  layers  appear 
in  localized  subdomains  of  fl  In  contrast,  for  problem  (1),  boundary  layers  ap¬ 
pear  throughout  the  perforated  domain  f2e.  Again  the  effect  of  the  convection  is 
expected  to  be  more  critical. 

Examples  of  applications  of  (1)  include  the  (Navier-)  Stokes  equation  in  per¬ 
forated  media.  The  twofold  difficulty  of  such  problems  lies  in  the  presence  of  two 
small  scales  that  we  detail  below.  Efficient  numerical  approaches  have  been  pro¬ 
posed  to  address  either  of  the  two  difficulties,  but  not  the  two  of  them  together: 

(i)  in  the  absence  of  convection  (i.e.  if  b  =  0  in  (1)),  the  problem  becomes 
a  purely  diffusive  problem  set  on  a  multiscale  perforated  domain.  This 
problem  can  be  efficiently  approximated  by  dedicated  MsFEM  approaches, 
such  as  the  one  we  proposed  in  [EOARD2]  (see  also  [EOARD5]),  which  is 
based  on  using  Crouzeix-Raviart  type  basis  functions.  See  Figure  2  for  a 
representative  example. 

(ii)  in  the  absence  of  perforations,  the  problem  reads 

—aAu  +  b{x )  ■  Vw  —  f  in  fl,  u  =  0  on  dQ,  (5) 

and  is  again  simple  to  treat  using  the  now  classical  stabilization  techniques 
for  finite  element  methods  (FEM).  For  example,  the  exact  solution  of  (5),  in 
the  one-dimensional  case,  is  shown  on  Figure  3  for  b  large.  It  is  well-known 
that,  unless  the  mesh  size  is  smaller  than  the  boundary  layer  width,  classical 
PI  FEM  perform  poorly.  The  solution  is  not  well  approximated,  not  only 
within  the  boundary  layer,  but  also  away  from  it,  because  of  numerical 
instabilities.  Methods  alternative  to  classical  PI  FEM,  known  as  stabilized 
approaches  (e.g.  the  Streamline  Upwind  Petrov- Galerkin  (SUPG)  method), 
have  been  introduced  long  ago  in  the  literature.  They  are  very  effective,  as 
can  be  seen  on  Figure  3. 
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Figure  2:  Left:  an  example  of  a  domain  with  non-periodic  perforations  (represented 
in  black).  Right:  relative  errors  with  various  approaches:  FEM  -  the  standard  Ql 
finite  elements,  no  OS  -  MsFEM  with  linear  boundary  conditions,  osc  -  MsFEM 
with  oscillatory  boundary  conditions,  OS  -  MsFEM  with  oversampling,  CR  -  the 
Crouzeix-Raviart  type  MsFEM  approach  we  propose.  Results  for  all  these  methods  are 
represented  by  solid  lines.  The  dashed  lines  correspond  to  the  variants  of  these  methods 
where  we  enrich  the  finite  element  spaces  using  bubble  functions.  See  [EOARD2]  for 
details. 


Figure  3:  In  blue,  the  exact  solution  to  (5)  for  17  =  (0, 1)  and  b  =  200.  In  green  (resp. 
red),  the  numerical  solution  obtained  by  a  PI  (resp.  stabilized  PI)  approach  with 
H  =  0.08.  The  PI  solution  shows  spurious  oscillations  while  the  PI  SUPG  solution  is 
very  accurate  except  in  the  boundary  layer. 


10 


DISTRIBUTION  A.  Approved  for  public  release:  distribution  unlimited. 


In  the  sequel,  we  focus  on  problem  (1).  The  question  is  to  design  an  efficient 
MsFEM  approach  for  that  problem.  It  is  unclear  whether  it  needs  to  be  stabilized 
and  how  to  achieve  this  in  the  most  effective  way.  We  report  below  on  the 
preliminary  results  that  we  have  obtained,  and  refer  to  [EOARD4,  EOARD7]  for 
comprehensive  results. 

1.2  Numerical  approaches 

We  recall  that  the  variational  formulation  associated  to  (1)  consists  in  finding 
u£  G  Hq(Q.s)  such  that 


Vt>  e  Ac(u‘,v)  =  F,(v), 


(6) 


where 


Ae(u,  v)  —  a  f  Vu  •  V'u  +  (bMA  .  Vu  j  v  and  Fe{y)  —  I  fv. 

Jne  \  £  J  Jne 

Following  the  spirit  of  the  MsFEM  approaches,  our  aim  is  to  introduce,  at  the 
coarse  scale,  a  finite-dimensional  approximation  space  spanned  by  suitably  cho¬ 
sen  basis  functions,  and  to  consider  the  Galerkin  approximation  of  (6)  in  that 
space.  We  describe  here  the  different  variants  we  have  considered,  before  turning 
in  Section  1.3  to  the  obtained  numerical  results. 

We  mesh  12,  and  denote  Th  the  set  of  triangles  K  C  12.  We  denote  by  Sh  the 
set  of  all  the  internal  edges  of  Th-  Note  that  we  mesh  12  and  not  the  perforated 
domain  12e.  This  allows  us  to  use  coarse  elements  (independently  of  the  hne  scale 
present  in  the  geometry  of  12e),  and  leaves  us  with  a  lot  of  flexibility.  The  mesh 
does  not  have  to  be  consistent  with  the  perforations  Bs.  Some  nodes  may  lie  in 
B£,  and,  likewise,  some  edges  may  intersect  Be.  We  assume  in  the  sequel  that  no 
element  K  e  Th  and  no  internal  edge  E  e  Sh  is  a  subset  of  the  perforations  Be 
(recall  that  the  mesh  is  coarse  while  the  perforations  are  small,  so  this  assumption 
is  easily  satisfied  in  practice). 

In  the  vein  of  our  earlier  work  [EOARD2]  on  MsFEM  approaches  for  diffusive 
problems  set  on  perforated  domains,  we  consider  MsFEM  basis  functions  defined 
with  Crouzeix-Raviart  type  boundary  conditions.  For  purely  diffusive  problems, 
they  have  proved  to  be  extremely  flexible. 

We  now  briefly  present  our  four  numerical  approaches.  For  each  of  them,  the 
basis  set  consists  of  two  types  of  functions: 
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•  a  function  e  £n  associated  to  each  edge  E  of  the  mesh,  as  is  usual  for 
Crouzeix-Raviart  type  approaches; 

•  a  function  vhe,K  associated  to  each  element  K  6  77/  of  the  mesh,  which 
vanishes  (in  a  weak  way)  on  the  edges  of  that  element.  Such  functions 
are  called  bubble  functions  in  what  follows.  In  the  context  of  purely  dif¬ 
fusive  problems,  the  need  to  use  such  basis  functions  for  problems  posed 
on  perforated  domains  has  been  detailed  in  our  previous  work  [EOARD2], 
and  can  be  observed  on  Figure  2.  In  the  present  context  of  convection- 
diffusion  problems,  we  have  also  considered  the  variant  where  we  do  not 
include  bubble  functions  in  the  basis  sets.  Results  are  shown  on  Figure  4  for 
the  Adv-MsFEM  variant  described  in  Section  1.2.1  below  (similar  results 
have  been  obtained  with  the  other  MsFEM  variants).  As  for  our  previous 
work  [EOARD2],  the  interest  of  adding  bubble  functions  to  the  basis  set  is 
evident. 


Figure  4:  Relative  H 1  error  (as  defined  by  (20)  below)  for  the  Adv-MsFEM  vari¬ 
ant  (a  =  1/4,  £  =  0.03;  the  vertical  dashed  line  represents  the  region  where 
H  =  e),  when  including  or  not  bubble  functions  in  the  basis  set. 
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1.2.1  Adv-MsFEM  variant 

In  that  variant,  for  each  K  e  Th,  we  solve  the  local  problem 

-«Ar’K  +  bS*M.  .  v^’K  =  1  ill  K  \  ~Bei  =  0  in  Be,  (7) 

£ 

with  the  so-called  Crouzeix-Raviart  boundary  conditions:  for  each  edge  Tj  of  K, 
we  set  /  TfK  =  0  and  n  ■  Vvh£,K  =  Xt  on  R  for  some  constant  A*.  The  basis 

J  it 

functions  are  called  bubbles  in  what  follows. 


Second,  for  each  internal  edge  E  E  £h,  we  denote  by  Kg  and  K|  the  two 
triangles  sharing  this  edge,  set  Kg  :=  (see  Figure  5),  and  consider  $>£,E 

solution  to 


-aA$£'E  +  =  0  in  Kg  \  Be, 

-aA$£,E  +  bS*Ill  .  =  0  in  K|  \  ~Be, 

£ 

=  0  in  B£, 


(8) 


with,  for  each  edge  E'  C  <9K e,  /  ®£,E  =  0  and  n  •  V$£’E  =  A e>  on  E'  for  some 

J  E' 


constant  A e>  and  /  <^>e’  =  1  and  n  ■  W^£’E  —  Xe  on  E  for  some  constant  A 

Je 

(with  an  a  priori  different  constant  on  the  two  sides  of  E). 


K 


2 

E 


Figure  5:  Construction  of  &£,E  following  (8). 


The  Adv-MsFEM  approximation  space  is 

V%Adv  =  Span  {$£’E,  T£’k,  E  E  Sh,  K  eTh}  ■  (9) 

The  Adv-MsFEM  method  is  the  standard  Galerkin  approximation  of  (6)  on 

t  7-£,Adv 

VH  : 

Find  u%  E  V%Adv  such  that,  \/v%  E  V%Adv ,  AE{u%,v %)  =  Fe(v %), 
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where 


AE(u,v)  =  a  ^  f  Vu  •  Vm  +  ^2  [  .  \7u  \  v_ 

K.£7jj  Knf2e  KgTh  V  e  / 

1.2.2  Classical  MsFEM  variant 

The  basis  functions  of  the  Adv-MsFEM  variant  presented  above  depend  on  the 
convection  held  b.  Thus,  if  b  changes,  the  basis  set  has  to  be  recomputed.  This 
may  be  expensive.  We  explore  here  the  possibility  of  a  MsFEM  variant  in  which 
the  basis  functions  do  not  depend  on  b.  This  variant  is  henceforth  called  the 
MsFEM  variant. 

For  each  K  e  Eh,  we  solve  the  local  problem 

-«ATo’k  =  1  in  K\  ~Be,  Tq’k  =  0  in  B£,  (10) 

with  the  same  boundary  conditions  as  for  (7).  In  addition,  for  each  internal  edge 
E,  we  consider  solution  to  (see  Figure  5  for  notations) 

f  -«A$g£  =  0  in  Kg  \ 

i-aA^  =  0inK|\Se,  (11) 

[%’E  =  0  in  B£, 

with  the  same  boundary  conditions  as  for  (8). 

The  MsFEM  approximation  space  is 

Vff  =  Span  To’K,  Ee£H,  K  e  TH)  •  (12) 

The  MsFEM  method  is  the  standard  Galerkin  approximation  of  (6)  on  Vfj\ 

Find  u£h  e  Vh  such  that,  \/veH  G  Vfj,  AE(u£H,v£H )  =  F£(y£H). 

In  the  numerical  results  described  below  (see  Section  1.3),  we  will  observe  that 
this  approach  performs  poorly  when  a  decreases,  that  is  when  convection  in¬ 
creasingly  dominates.  There  is  thus  a  need  for  stabilization. 
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1.2.3  Stab-MsFEM  variant 

We  use  the  same  approximation  space  VR  as  in  the  classical  MsFEM  variant 
of  Section  1.2.2,  but  we  modify  the  bilinear  form  by  adding  stabilization  terms. 
The  problem  consists  in  finding  ueH  G  V R  such  that 

VyH  £  Vh h  (u£h,v £h)  +  Astab(u£H , v£h)  =  F£{v£H )  +  Fsta b(vjj),  (13) 


where 


A-stab  (u%  ,v%)  — 

K gTh  ' 

/  tk(x)  Cu£h 

/Knoe  L  J  . 

bix/e)  _  p 
•  VveH 

£ 

(14) 

-^stab  (' VH )  —  y  ' 

KP  TV  ' 

/  7*  (X)fb(XK 

/Knne  £ 

<1 

c* 

(15) 

In  the  definition  of  *4.stab  and  Fstab,  rK  is  a  stabilization  parameter,  the  value  of 
which  is  here  chosen  to  be 

TkW  =  2\b(xH/e)/e\  (COth(Pe(l) H)  ~  '  Pe(x)  = 

and 

„  .  b(x/e)  „ 

Cu  =  —aAu  H - -  V« 

£ 

is  the  operator  corresponding  to  the  problem  at  hand. 

The  method  is,  as  is  well  known,  strongly  consistent.  In  view  of  the  defini¬ 
tion  (10)— (11)— (12)  of  VR,  any  ueH  G  VR  can  be  written  as 

<4M=  T  ue%-b(x)+  V  c/k^k(i), 

E€l£h 

and  thus,  for  any  ueH  and  veH  in  VR,  we  have  that 

•Astab  (^if  iV£jj)  =  Aup  w{Wjj,Vjj)  (16) 

where 
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In  practice  however,  we  only  know  discrete  approximations  'kg  and  d’gf ,  on  a 

fine  mesh  K/,,  of  the  solutions  vIfij'K  and  $0^  to  (10)  and  (11),  respectively.  Put 
differently,  we  manipulate 

VS,„  =  Span  {»*£ ,  9'*,  Ee£„,  K  €  T„  } 


instead  of  Vfj.  It  follows  that,  when  we  use  a  P1  approximation  on  a  fine  mesh 
K/j  for  the  local  problems  (10)  and  (11),  may  be  discontinuous  at  the 

edges  of  the  mesh  K/-(,  thus  —A u£Hh  ^  ^^(K)  and  eventually  Astabiu£Hh)  veH  h) 
is  ill-defined. 

We  may  consider  at  least  two  ways  to  circumvent  that  difficulty.  First,  we 
may  define  the  stabilization  term  as 


A. 


stab 


,uH,hi  VH,h ) 


E  E 

KeTk  KCKft 


/  TK(x) 

Mkjh 

//■vPlQe 

l  J 

_ 

b[xje ) 


We  then  obtain  a  strongly  consistent  stabilized  method.  We  will  however  not 
proceed  in  this  direction  and  favor  an  alternate  approach,  to  which  we  now  turn. 

Based  upon  the  observation  (16)  for  the  ’’ideal”  space  V§,  we  may  use  the 
stabilization  term  (17)  rather  than  (14).  In  contrast  to  (14),  the  quantity  (17)  is 
also  well  defined  on  V^h.  The  Stab-MsFEM  variant  we  employ  is  hence  defined 
by  the  following  variational  formulation:  find  ueH  h  G  h  such  that 


e  Y£|fc, 


{u£Hh,  v£H  h)  +  Anpw(u£H  h,  v£H  h)  —  Fe{v£H  h )  +  Fs tab{v£H  h).  (18) 


We  emphasize  that  employing  that  stabilization  comes  at  a  price:  we  give  up  on 
strong  consistency.  We  note  that  we  have  already  proceeded  likewise  in  [EOARD3], 
and  that  we  were  able  to  prove  there,  in  some  cases,  that  the  method  is  conver¬ 
gent  despite  the  absence  of  consistency. 


1.2.4  Stab-Adv-MsFEM  variant 

In  echo  to  the  variant  described  in  Section  1.2.3,  it  may  be  interesting  to  also  con¬ 
sider  a  stabilized  version  of  the  Adv-MsFEM  variant  introduced  in  Section  1.2.1. 
The  problem  consists  in  finding  ueH  G  V^Adv  such  that 

Vv£h  G  V^A%  A£(u£h,  v£h )  +  Atab(?4,  v£H )  =  F£{v£h)  +  Fstab(A),  (19) 

where  V^Adv  is  defined  by  (9)  and  Atab  and  Fstab  are  defined  as  in  (14)-(15). 
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In  practice,  we  again  only  know  discrete  approximations  \I/^K  and  ,  on 
a  fine  mesh  K/t,  of  the  solutions  Te’K  and  Qe,E  to  (7)  and  (8),  respectively. 
The  term  Astab(u£H  h,  veHh)  is  again  ill-defined.  We  proceed  as  for  the  Stab- 
MsFEM  variant  (see  (18))  and  consider  the  following  variational  formulation: 
find  u£Hh  e  V^dv  such  that 


VvH,h  e  ^H,h  )  Af  (u£H  h,  v£H  h)  +  Anpw(u£Hh,  v£H  h)  —  Fe(v£H  h)  +  Fstah(v£Hh), 
where 

AnPw(ueH}h,v£Hh)  =  Y]  UK  /  tk(x) 

./KnOe 


K  eTH 


b(x/e) 


■  Vv£ 


with  u£Hh{x)  —  y  UE$hE(x)+  y  UKFehK(x). 

EeSjj  KgT h 


1.3  Numerical  results 

We  describe  here  some  preliminary  numerical  tests,  and  refer  to  [EOARD4]  for 
more  comprehensive  results. 

We  consider  problem  (1)  in  dimension  d  —  2.  The  perforated  domain  is 

=  12  \  B£,  where  12  =  (0,  l)2  and  B£  is  the  set  of  squares  of  sidelength  0.5£ 
periodically  located  on  the  regular  grid  of  period  e.  We  consider  the  generic 
case  when  the  mesh  edges  intersect  the  perforations.  Unless  otherwise  stated, 
we  set  b  =  (1, 1)T  and  /  =  1.  The  parameters  £  and  a  will  be  specified  along  the 
different  numerical  results. 

We  focus  on  the  relative  error  in  the  Ft1  norm,  which  is  defined,  for  any 
numerical  solution  v,  by 


elv  = 


E 


KeTff 


\v  Mrefllr/i(K\B7) 


(20) 


where  ure f  is  the  reference  solution,  which  has  been  computed  by  directly  solv¬ 
ing  (1)  on  a  very  fine  mesh. 


1.3.1  Comparison  of  the  MsFEM  approaches 

We  compare  the  four  variants  we  have  introduced,  namely  the  MsFEM  (see 
Section  1.2.2),  the  Stab-MsFEM  (see  Section  1.2.3),  the  Adv-MsFEM  (see  Sec¬ 
tion  1.2.1)  and  the  Stab-Adv-MsFEM  (see  Section  1.2.4).  All  the  methods  share 
comparable  offline  costs,  and  identical  online  cots. 
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We  first  investigate  the  sensitivity  to  the  coarse  mesh  size  H .  We  thus  fix 
a  =  1/4,  e  =  1/0.03  and  let  H  vary.  We  plot  the  errors  on  Figure  6.  We  observe 
that  the  stabilization  is  not  useful  (in  the  sense  that  the  methods  with  or  without 
the  stabilization  terms  yield  similar  results),  and  that  the  Adv-MsFEM  method 
provides  more  accurate  results  than  the  MsFEM  method. 


H 


MsFEM 
Stab-MsFEM 
♦  Adv-MsFEM 
Stab-Adv- MsFEM 


Figure  6:  Relative  error  (20)  for  the  4  numerical  variants  (a  —  1/4,  e  —  0.03;  the 
vertical  dashed  line  represents  the  region  where  H  =  e). 


We  next  investigate  the  influence  of  the  parameter  a,  which  measures  the 
relative  importance  of  the  diffusion  and  the  convection  in  the  problem.  We  fix 
H  =  1/16,  e  =  1/32  and  let  a  in  (1)  vary.  We  plot  the  errors  on  Figure  7.  Again, 
the  stabilization  does  not  appear  to  be  useful,  and  the  Adv-MsFEM  method 
provides  more  accurate  results  than  the  MsFEM  method. 


We  eventually  examine  the  robustness  of  the  approaches  to  the  small  scale 

£  present  in  the  problem.  MsFEM  approaches  are  indeed  supposed  to  provide 

accurate  results  for  fixed  H  even  if  £  is  asymptotically  small,  in  contrast  to 

/7ix\  fny\ 

standard  FEM  approaches.  We  set  f(x,y )  =  sin  sin  ^  = 

a  =  1/16  and  let  £  vary.  As  for  the  previous  tests,  we  observe  on  Figure  8  that 
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Figure  7:  Relative  error  (20)  for  the  4  numerical  variants  (e  =  1/32,  H  =  1/16). 


the  Adv-MsFEM  and  the  Stab-Adv-MsFEM  methods  provide  more  accurate 
results  than  the  MsFEM  and  the  Stab-MsFEM  methods. 

1.4  Conclusions  and  perspectives 

Our  first  numerical  tests  show  that  the  Adv-MsFEM  and  the  Stab-Adv-MsFEM 
approaches  provide  solutions  with  the  best  accuracy,  in  term  of  H 1  error  over  the 
whole  domain  However,  our  results  are  only  preliminary,  and  much  remains 
to  be  understood. 

For  instance,  to  have  a  better  understanding  of  the  different  numerical  ap¬ 
proaches,  we  wish  to  investigate  the  question  of  where  in  the  error  is  essen¬ 
tially  located.  This  is  a  delicate  question.  For  MsFEM  approaches,  the  error 
typically  comes  from  a  mismatch  between  the  exact  solution  and  the  numeri¬ 
cal  solution  close  to  the  edges  of  the  coarse  mesh  (thus  the  introduction  of  the 
popular  oversampling  variant,  ...).  However,  for  convection-dominated  prob¬ 
lems,  the  situation  can  be  different.  For  instance,  we  considered  in  [EOARD3] 
a  multiscale  convection-dominated  problem,  for  which  all  numerical  approaches 
produced  large  errors  (overshadowing  the  error  at  the  coarse  mesh  edges)  in  a 
boundary  layer  due  to  the  strong  convection.  For  the  problem  considered  here, 
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Figure  8:  Relative  error  (20)  for  the  4  numerical  variants  (a  =  1/16,  H  =  1/16). 


which  is  convection-dominated  and  set  on  a  perforated  domain,  the  situation 
might  be  again  different.  This  question  thus  has  to  be  investigated  in  details. 

Another  question  we  wish  to  study  is  to  better  understand  why  the  Adv- 
MsFEM  approach  appears  here  to  be  the  best,  whereas  it  used  to  be  outper¬ 
formed  by  the  Stab-MsFEM  approach  in  the  work  [EOARD3].  Note  that  criteria 
of  performance  are  different  in  both  cases:  the  H 1  error  in  the  former  case,  the  H 1 
error  outside  a  boundary  layer  in  the  latter  case  (due  to  the  fact  that,  as  pointed 
out  above,  all  numerical  approaches  produced  large  errors  in  the  boundary  layer). 

ft  also  remains  to  understand  whether  the  Stab-Adv-MsFEM  approach  pro¬ 
vides  systematically  better  results  than  the  Adv-MsFEM  approach.  If  not,  then 
the  Adv-MsFEM  approach  appears  to  be  the  method  of  choice,  since  it  provides 
accurate  results  and  does  not  require  the  assembling  of  stabilization  terms. 

All  these  questions  will  be  investigated  in  [EOARD4], 


20 


DISTRIBUTION  A.  Approved  for  public  release:  distribution  unlimited. 


2  Construction  of  coarse  approximations  for  prob¬ 
lems  with  highly  oscillatory  coefficients 

The  work  we  now  describe  addresses  the  construction  of  coarse  approximations 
for  problems  with  highly  oscillatory  coefficients.  It  has  been  completed  in  collab¬ 
oration  with  Simon  Lemaire,  a  postdoc  student  fully  funded  by  this  Grant  from 
June  1st,  2014  until  September  30th,  2015. 

Consider  an  heterogeneous  material,  modelled  by  highly  oscillatory  coeffi¬ 
cients,  which  we  denote  here  by  A£ .  For  simplicity,  we  assume  that  Ae  is  a 
symmetric  matrix.  The  question  we  investigate  is  to  find  the  best  non-oscillating 
coefficient  A  (think  e.g.  of  a  constant  coefficient)  that  is  consistent,  in  a  sense  to 
be  made  precise,  with  the  behavior  of  the  heterogeneous  material  modelled  by 
Ae .  Put  differently,  we  wish  to  approximate  the  solution  u£  to  the  problem 

—  div  [H£(a;)V«e]  =  /  in  Q,  u£  =  0  on  dfl,  (21) 

using  a  similar  problem  with  a  constant  coefficient  A,  which  reads 

—  div  [AVm]  —  f  in  ft,  u  —  0  on  <9f2.  (22) 

We  want  the  approximation  to  be  accurate  for  ideally  all  right  hand  sides  /. 

The  motivation  for  this  approach  is  to  find  a  pathway  alternative  to  stan¬ 
dard  homogenization  techniques  when  these  latter  are  difficult  to  use  in  practice. 
Indeed,  homogenization  theory  (and  the  numerical  homogenization  techniques, 
such  as  MsFEM,  that  derive  from  it)  requires  a  complete  knowledge  of  the  oscil¬ 
latory  coefficient  A£  (e.g.  to  compute  the  multiscale  basis  functions  in  (7)-(8)). 

In  many  practical  cases,  this  coefficient  is  only  partially  known,  or  one  only  has 
access  to  the  solution  of  (21)  for  some  loadings  /. 

When  £  is  asymptotically  small,  A  should  then  be  a  good  approximation  of 
the  homogenized  matrix  A*.  This  is  indeed  the  case  for  the  two  formulations, 
see  (23),  (24)  and  (25)  below,  that  we  have  considered  (see  [5,  EOARD1]).  The 
approach  can  thus  be  employed  to  approximate  A*.  However,  the  most  inter¬ 
esting  regime  is  when  e  is  small  but  not  asymptotically  small.  Our  approach 
provides  a  way  to  approximate  the  solution  u£(f)  to  (21)  (which  is  expensive 
to  compute,  especially  in  a  multi-query  context  where  several  /  have  to  be  con¬ 
sidered)  by  the  solution  u(f)  to  (22),  which  can  be  computed  in  an  inexpensive 
way  as  soon  as  A  has  been  identified.  Besides  A*  or  u£(f),  another  quantity  of 
interest  is  V«£(/).  We  examine  below  how  to  approximate  these  quantities. 
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2.1  A  quadratic  optimization  formulation 

A  central  question  is  to  construct  in  an  efficient  manner  the  best  constant  coeffi¬ 
cient  A  approximating  the  oscillatory  problem.  We  have  addressed  this  question 
in  [5],  where  we  have  considered  the  following  optimization  problem: 

inf  _  sup  Je(A,f)  (23) 

symmetric  matrix  A  f£L2(fl),  ||/||z,2(n)  =  l 

with 

Je(A,f)  =  \\u£(f)  -u(AJ)\\2L2{n)  (24) 

where  u£(f)  is  the  solution  to  (21)  with  the  right-hand  side  /  and  u(A,  f )  is  the 
solution  to  (22)  with  the  right-hand  side  /  and  the  constant  matrix  A.  However, 
this  problem  is  difficult  to  solve,  because  it  is  nonconvex.  As  pointed  out  in  the 
previous  EOARD  report  [EOARD6],  we  have  recently  investigated  a  different 
formulation  of  the  problem.  It  is  given  by  (23),  where  now  we  take 

MAf)  =  ||(-A)-1(divpV«'(/))  +  /)|[s(n),  (25) 

where  ue(f)  is  the  unique  solution  to  (21)  and  the  operator  (— A)-1  is  defined 
by 

w  =  (—A)~1g  when  —  Aw  =  g  in  12  and  w  =  0  on  <912. 

Note  that  the  cost  function  (25)  is  related  to  (24)  through  the  application,  inside 
the  L2  norm,  of  the  operator  (— A)-1  ^  div(AV-)^ . 

We  emphasize  that  (25)  is  quadratic  with  respect  to  /  and  A,  in  contrast 
to  (24).  The  problem  (23)-(25)  thus  has  the  advantage  of  being  convex.  This  is 
of  paramount  importance  for  the  efficiency  of  the  numerical  approach  [EOARD  1], 

2.1.1  Algorithm  to  solve  (23)— (25) 

An  algorithm  to  solve  (23)-(25)  has  been  presented  in  the  previous  EOARD 
report  [EOARD6].  We  simply  recall  here  that  the  supremum  over  /  G  L2(12) 
in  (23)  is  approximated  by  a  supremum  in  a  finite-dimensional  space.  More 
precisely,  problem  (23)  is  approximated  by 

inf  _  sup  Je(A,f ),  (26) 

symmetric  matrix  A  feVp,  ||/||i2(n)=1 
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where  Vp  is  the  finite  dimensional  space  generated  by  the  first  P  eigenfunctions 
of  the  laplacian  operator  on  fl,  that  is 

VP  =  Span  {fp,  1  <  p  <  P}  ,  (27) 

where  —A fp  =  A pfp  on  with  fp  —  0  on  <9fl  We  normalize  the  functions  fp  so 

that  I  fpfq  &qp- 

Jn 


2.1.2  Adaptation  to  the  random  setting 


We  have  presented  in  the  previous  EOARD  report  [EOARD6]  some  preliminary 
numerical  results  in  a  periodic  case.  In  this  report,  we  focus  on  the  stationary  ran¬ 
dom  case.  Our  approach  is  adapted  to  that  context  by  simply  considering  (26), 
where  Je(A,  f)  is  now  given  by 


UAf) 


(— A)-‘(div  (A VE [«'(/)])+/) 


2 

L2(Q)  ’ 


where  u£(uj,  f )  is  the  unique  solution  to 


(28) 


div 


A(^,w)  Vu£ 


=  /  in  0,  u£  =  0  on  d 0, 


(29) 


where  the  matrix  A  is  random  and  stationary  (i.e.  statistically  homogeneous). 

In  practice,  the  expectation  in  (28)  is  approximated  by  an  empirical  mean 
over  M  i.i.d.  realizations: 


M 


E  [«*(/)] 


M 


'y  ^  u  iyimi  /)• 


m=  1 


Problem  (26)  is  thus  approximated  by 

( 


< 


inf  _  sup  J£  (A ,/), 

symmetric  matrix  A  f£VP,  ||/||£2(fJ)=l 

M 


-A)-‘  |  div  [  Av 


y  ^ u  iyjm  i  f ) 

m=l 


+  / 


(30) 


L2(Q) 


In  the  sequel,  we  compare  our  approach  with  the  standard  homogenization 
approach,  which  reads  as  follows.  The  solution  u£(-,uj )  to  (29)  converges,  when 
£  — y  0,  to  M* ,  solution  to  the  homogenized  problem 

—  div  [A*V«*]  =  /  in  Q,  u*  =  0  on  d£l,  (31) 
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where  the  homogenized  matrix  A*  is  constant  and  deterministic.  This  matrix  is 
challenging  to  compute,  and  in  practice  approximated  as  follows.  Consider  the 
truncated  corrector  equation 

—  div  (A(-,u)(p  +  Vwp  (-,u;)))  =  0  in  QN ,  (-,u)  is  Qh -periodic,  (32) 

posed  on  a  large  domain  QN  =  (-N,  N)d ,  where  p  6  Rd  is  a  constant.  We  then 
compute  the  random  matrix  A*(uj)  defined  by 


KM] 


(33) 


where  {ej}1<i<d  is  the  canonical  basis  of  M.d.  It  is  well-known  (see  e.g.  [4])  that 
A*  (ta)  almost  surely  converges,  when  N  — >  +oo,  to  A*.  Since  A*  (uj)  is  random, 
it  is  natural  to  consider  M  independent  and  identically  distributed  (i.i.d.)  re¬ 
alizations  of  the  held  A,  say  {A (*,  wm)}1<m<M,  solve  (32)  and  compute  (33)  for 
each  of  them,  and  define 


an,m 


m=  1 


(34) 


as  a  practical  approximation  to  A*. 


2.1.3  Approximation  of  u£  in  the  U1-norm 

We  have  pointed  out  above  that  our  approach  yields  a  matrix  Ae  which  is  a 
converging  approximation  of  A*.  The  solution  to  (22)  with  A  =  Ae,  which  we 
denote  u(Ae),  is  thus  a  converging  approximation  (in  the  H 1  norm,  and  thus  in 
the  L 2  norm)  of  w*.  Since  u£(-,  u)  converges  (in  the  L 2  norm)  to  «*,  we  therefore 
see  that  our  approach  allows  us  to  obtain  a  converging  approximation  of  u£(-,(d). 
However,  in  the  H 1  norm,  u£(-,u)  is  not  close  to  «*,  and  hence  u£(-,co)  and  u(Ae) 
are  not  close  to  each  other. 

We  explain  here  how  our  approach  can  be  complemented  to  obtain  an  ap¬ 
proximation  of  E  [w£]  in  the  171-norm. 

In  many  settings  of  homogenization  theory  (and  in  particular  in  the  ran¬ 
dom  setting  we  consider  here),  once  the  corrector  problems  (in  practice,  prob¬ 
lems  (32))  are  solved  to  compute  the  homogenized  matrix,  the  first-order  two- 
scale  expansion  approximates  u£  in  the  i/1-norm.  Consider  indeed  the  function 

U£'\x,  u)  =  U*(x)  +  £  ^2  Wei  (*,w) 

i= 1  ~  1 
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where  wei  is  the  unique  solution  with  vanishing  mean  value  to  the  corrector 
equation  for  p  =  e*.  It  is  known  that,  in  average,  u£  —  u£}1  converges  to  0  in  the 
H1  norm,  in  the  sense  that  (see  [6,  Theorem  3]) 


lim  E 

£ — ^0 


u 


S^) 


u 


£,1, 


■,uj) 


2 


0. 


This  implies  that 


E  [V«£(-,  a;)]  =  C£  V«*  +  h.o.t., 


(35) 


where  Ce  is  a  d  x  d  oscillatory  matrix  field  given  by,  for  any  1  <  i,  j  <  d, 


[C'Ji.i  =  1  +  E  [diWei(-/£,u)\, 


[Ce]itj  =  E  [di.wej  (■/£,  w)]  if  j  ^  i.  (36) 


Onr  idea  for  constructing  an  approximation  of  E  [Vu£]  is  to  mimick  for¬ 
mula  (35)  and  seek  an  approximation  under  the  form  Ce'VTie,  where  u£  =  u(Ae). 
Once  the  best  matrix  Ae  has  been  computed,  we  compute  a  surrogate  C£  of  C£ 
by  solving  the  least-squares  problem 


inf 

Ce(L2(Q))dxd  ^ 


VEWH1-C  Vue(/r)||L(n)'‘ 


(37) 


for  a  given  number  R  of  right-hand  sides. 

In  practice,  the  expectation  in  (37)  is  approximated  by  a  mean  over  M  re¬ 
alizations,  and  the  right-hand  sides  fr  selected  for  (37)  are  the  first  R  basis 
functions  of  the  space  Vp  defined  by  (27),  with  R  such  that 


R<  P. 


The  fact  that  we  consider  a  number  R  <  P  of  right-hand  sides  makes  the  H1- 
reconstruction  an  inexpensive  post-processing  procedure  once  the  best  matrix  is 
computed,  as  we  already  have  at  our  disposal  ue(oj,  fr)  for  1  <  r  <  R. 

The  tensor  C  is  discretized  as  a  piecewise  constant  function  on  a  fine  mesh. 
Solving  (37)  hence  amounts  to  solving  many  least-squares  problems  (independent 
one  from  each  other)  set  in  a  d  x  d  dimensional  space  (as  C  is  a  d  x  d  matrix). 
The  associated  cost  is  negligible  in  comparison  to  the  cost  for  computing  the 
solutions  to  (29)  that  are  needed  in  (30). 
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2.2  Numerical  results 


We  work  in  dimension  d  —  2,  with  =  (0,  l)2.  We  consider  the  random  checker¬ 
board  test-case,  namely 


Ae(x,  y,  to)  =  asta(x/£,  y/e}  u)  Id2,  (38) 

with  asta  given  by 

asta{x,  y,  u)  =  Y  1Q+k(x,y)Xk(u),  (39) 

fcez2 

where  Q  =  (0,  l)2  and  where  the  random  variables  X &  are  i.i.d.  and  such  that 
P(A"fc  =  4)  =  P(A"fc  =  16)  =  1/2.  An  explicit  expression  for  the  homogenized 
matrix  is  known  in  that  case  (of  course,  our  approach  does  not  use  that  knowl¬ 
edge): 

A*  =  8  Id2.  (40) 

In  what  follows,  we  only  consider  the  case  of  small  values  of  £  in  (29)  (namely, 
£  <  | Q | 1/2/ 10  =  0.1).  We  refer  to  [EOARD1]  for  comprehensive  results. 

In  short,  the  numerical  experiments  reported  below  show  that  the  approxi¬ 
mation  of  A*  obtained  by  the  classical  homogenization  approach  is  slightly  more 
accurate  than  the  one  obtained  with  our  approach.  In  contrast,  our  approach 
provides  a  better  L2- approximation  and  a  better  H 1  -approximation.  In  terms  of 
computational  cost,  our  approach  is  less  expensive  for  moderately  small  values 
of  e,  and  more  expensive  for  asymptotically  small  values  of  e. 

2.2.1  Cost  comparison  to  the  standard  homogenization  approach 

For  asymptotically  small  values  of  £,  our  method  can  be  seen  as  a  practical 
approach  to  compute  the  homogenized  matrix  A*.  The  question  is  whether 
this  approach  is  efficient,  in  comparison  with  the  classical  approach  in  random 
homogenization.  As  recalled  above,  this  latter  approach  requires  solving  the 
equations  (32),  set  on  QA!  =  (—N,N)d.  The  coefficients  of  (32)  vary  at  scale 
1.  In  that  case,  to  hope  for  an  accurate  approximation,  one  has  to  consider  a 
meshsize  1.  In  our  approach,  we  need  to  solve  (29),  where  the  coefficients 

vary  at  scale  e.  To  that  aim,  we  use  a  mesh  of  size  h.  We  choose  h  and  H  such 
that 

£  1 
h  =  J{' 
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hoping  thus  for  a  similar  accuracy  when  solving  (29)  or  (32). 

We  see  that,  up  to  an  appropriate  choice  of  the  parameter  H  such  that 


/2  N\d  size(fl) 
hd  ’ 


(42) 


both  the  standard  homogenization  approach  and  our  approach  require  solving 
linear  systems  of  the  same  size.  The  computational  workload  for  the  two  ap¬ 
proaches  is  thus  of  the  same  order  of  magnitude,  although  not  identical.  In  the 
sequel,  we  have  enforced  (42)  and  (41),  which  imply  that  N  in  (32)  and  e  in  (29) 
are  related  by 

N  =  size(fl)/2£.  (43) 

For  the  two  methods,  the  same  number  M  of  Monte  Carlo  realizations  is  used, 
and  the  same  M  realizations  are  considered. 

For  the  standard  approach  in  homogenization,  we  denote  by  A*’^1  the  practi¬ 
cal  approximation  of  A*)M  defined  in  (34).  Our  approach  consists  in  computing 
the  best  matrix  A£’h  following  (30).  Recall  that  P  denotes  the  dimension  of  the 
set  Vp  used  to  approximate  the  space  L 2  in  the  sup  problem  (see  (27)).  In  the 
sequel,  we  take  P  =  3. 


We  first  comment  on  the  CPU  times  (see  Figure  9),  namely  the  time  needed 
to  compute  either  A by  standard  random  homogenization  or  Ae  ’h  using  our 
approach.  To  perform  that  comparison,  we  make  use  of  an  implementation  that 
does  not  exploit  parallelism,  and  we  solve  the  linear  systems  by  means  of  an 
iterative  solver.  In  view  of  Figure  9,  our  method  is  slightly  faster  than  standard 
random  homogenization  for  values  of  N  up  to  approximately  14. 

This  observation  can  be  explained  as  follows.  For  the  number  M  =  100  of 
Monte  Carlo  realizations  that  we  consider,  we  can  neglect,  in  our  procedure, 
the  cost  of  optimization  in  comparison  to  the  cost  of  computing  the  solutions 
u£(um,fp )  to  (29),  for  1  <  m  <  M  and  1  <  p  <  P.  Hence,  to  compute  Ae'h  , 
we  have  to  (i)  assemble  M  =  100  stiffness  matrices,  (ii)  assemble  P  =  3  right- 
hand  sides,  and  (iii)  solve  P  x  M  —  300  linear  systems.  In  contrast,  to  compute 
A+'h,  one  has  to  solve  d  x  M  —  200  approximate  corrector  equations  (32),  that 
is  to  say  (i)  assemble  M  =  100  stiffness  matrices,  (ii)  assemble  d  x  M  —  200 
right-hand  sides,  and  (iii)  solve  d  x  M  —  200  linear  systems.  Consequently,  on 
the  one  hand,  our  approach  necessitates  solving  100  more  linear  systems  than 
the  standard  homogenization  approach,  but  on  the  other  hand,  the  standard 
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approach  necessitates  assembling  approximately  200  more  right-hand  sides  than 
our  approach.  This  explains  what  we  observe.  When  the  value  of  N  is  not  too 
large,  the  assembly  cost  is  higher  than  the  inversion  cost,  and  our  approach  is 
faster. 


Figure  9:  Ratio  between  the  CPU  time  needed  by  our  approach  and  the  CPU 
time  needed  by  the  standard  approach  in  random  homogenization  in  function  of 
N,  for  M  =  100. 


2.2.2  Approximation  of  A* 


We  consider  the  parameters  (W)3<fc<5  (defining  the  domain  on  which  we  solve  the 
corrector  problems  (32))  such  that  Nk  =  2k .  In  agreement  with  (43),  this  implies 
that  the  small  scale  present  in  (29)  is  given  by  ( £k)3<k<5  such  that  £k  =  2~(k+1\ 
The  associated  meshsizes  ( hk)s<k<5  and  ( Hk)z<k<b  are  computed  respectively 
letting  hk  =  e^/r  for  r  =  27,  r  =  54  or  r  =  108  and  using  (42).  We  consider 
M  =  100  Monte  Carlo  realizations. 


The  error  in  the  approximation  of  the  homogenized  matrix  is  defined,  for 

AM  e  {<#,  X:h),  by 


err  mat  = 


E 


AM 


h3 


~  [A*] 


hi 


2  \  1/2 


where  A*  is  taken  equal  to  the  exact  value  (40). 
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The  numerical  results  are  collected  in  Figure  10.  We  observe  that  the  ma- 
trix  Ae’h  provided  by  our  approach  converges  to  the  homogenized  matrix  as 
N  =  l/(2e)  grows.  However,  for  any  value  of  N  in  the  range  we  consider,  the 
approximation  of  A*  obtained  with  the  usual  approach  is  slightly  more  accurate 
(and  less  expensive  when  N  >  14,  see  Figure  9)  than  the  one  obtained  with  our 
approach. 


Figure  10:  Approximation  of  A*  by  standard  random  homogenization  (blue)  and 
by  our  approach  (black)  in  function  of  N,  for  M  =  100  realizations.  Since  M  is 
finite,  we  compute  the  error  err_mat  100  times.  The  thick  line  corresponds  to 
the  mean  value  over  the  100  computations  of  the  error.  The  dashed  lines  show 
the  95%  confidence  interval.  Results  obtained  with  h  such  that  s/h  =  27  (resp. 
e/h  =  108)  are  denoted  with  x  (resp.  o). 


2.2.3  Approximation  of  u£  in  the  L 2  norm 

We  denote  by 

•  ue!h(f)  empirical  expectation,  defined  by 

M 

U^h(f)  =  JJ  Ue,h(Um,  /),  (44) 

m=  1 

of  the  discrete  solutions  to  (29)  with  oscillatory  coefficients  given  by  (38)- 
(39)  and  right-hand  side  /; 
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•  u^hif)  the  discrete  solution  to  (31)  with  exact  matrix  (40)  and  right-hand 
side  /  (note  that  the  exact  homogenized  matrix  A*  is  usually  unknown); 

•  u^f  (f  )  the  discrete  solution  to  (31)  with  matrix  A*'^1  and  right-hand  side 


•  uf’^/(/)  the  discrete  solution  to  (22)  with  matrix  and  right-hand  side 
/•’ 

So  as  to  assess  the  quality  of  the  approximation  of  E[w£]  by 

uh{f)  e  {«*,/»(/),  <^(/),  «5f(/)} 

in  the  L2  norm,  we  define  the  criterion 


/ 


err  12  = 


sup  \\u£)h{f)  -  Uh(f) 
lev's 


1 2  \  V2 

lL2(n) 


u 


M 

e,h 


(fe 


2 

L2(T2) 


where  f£  E  Vq  denotes  the  argument  of  the  sup  problem  in  the  numerator.  Note 
that  the  supremum  is  taken  over  f  E  Vq  with  Q>P.  We  take  Q  =  16,  and  we 
have  checked  that  our  results  do  not  change  for  a  larger  value  of  Q. 


The  numerical  results  are  collected  in  Figure  11.  We  observe  that  the  solution 
<f(/)  provided  by  our  approach  is  a  better  L2-approximation  (for  the  range 
of  parameters  considered  here)  of  E(we)  than  the  solutions  associated  with  the 
exact  or  approximate  homogenized  matrices.  Due  to  the  small  number  of  right- 
hand  sides  we  consider  to  compute  Ae  ’h  ,  this  good  accuracy  is  not  an  immediate 
consequence  of  our  practical  procedure.  We  also  observe  that  the  accuracy  of  the 
three  approximations  and  u£'^'  improves  when  h  decreases  or  when 

M  increases,  in  somewhat  a  complex  manner.  In  terms  of  cost,  our  approach  is 
again  less  expensive  than  the  classical  approach  for  N  <  14. 


2.2.4  Approximation  of  u£  in  the  H 1  norm 


We  eventually  turn  to  the  W1-error.  We  denote  by  C(f^M  the  practical  approxi¬ 
mation  of  the  deterministic  oscillatory  matrix  Ce  defined  by  (36)  by  an  empirical 
mean  over  M  realizations  of  the  truncated  corrector  functions,  solutions  to  (32): 


nN,M 


ij 


1 


M 

3 ij  “h  ^  ^  /£,UJm). 

m=  1 
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Figure  11:  Approximation  of  E(m£)  in  the  L2-norm  (err_12  )  by  (red), 

(blue)  and  u X  (black)  in  function  of  N  (curves  with  x:  e/h  =  27  and  M  =  100; 
curves  with  o:  e/h  =  108  and  M  =  100;  curves  with  +:  e/h  =  27  and  M  =  400; 
curves  with  □:  e/h  =  54  and  M  =  400). 


For  /  G  Z/2 (0) ,  we  denote  by  and  C£f/u'Vu1/1//1  (f)  the  two  discrete 


N,Mi 


N,M  / 


e,h 


e,h 


equivalents  of  Ce  Vn*(/),  the  homogenization-based  approximation  of  E  (V«£(/)), 
obtained  by  using  the  exact  homogenized  matrix  (40)  and  the  matrix  re¬ 

spectively,  to  compute  an  approximation  of  «*(/).  In  our  approach,  we  seek 
a  discrete  approximation  of  E(Vu£(/))  under  the  form  C£’h  Vnff(/),  with 
R  =  P  =  3  (see  (37)).  For 


C“hVuh(f)  e  [C^M  V«*lfc(/),  CX1  (/),  CQT  V<f(/)j , 

we  define  the  criterion 


-yN,M  N,M  /  r\  ~^R,M  ^7 — P,M  / 


err  hi  = 


sup 

/evs 


V<fc(/)  -  V«fc(/) 


M 


\ 


L2(0) 


V«"fc(/e) 


L2(n) 


(45) 


where,  here  again,  the  supremum  is  taken  over  the  space  Vq  for  Q  =  16  P 
and  where  fe  G  Vq  denotes  the  argument  of  the  sup  problem.  We  recall  that, 
in  (45),  uXf)  is  the  empirical  mean  (44)  over  M  realizations  of  uejh{f',w).  It  is 
thus  an  approximation  to  E  [u£(f)]. 
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The  numerical  results  are  collected  in  Table  1.  Our  surrogate  defines  an 
approximation  of  E (Vu£)  which  is  systematically  more  accurate  than  the  one 
provided  by  the  homogenization-based  approach,  for  any  h  and  M. 


N 

8 

16 

32 

err_hl  for  C‘///Vl  'Vu+j,  ( e/h  =  27,  M  =  100) 

1.043  10"1 

(e/h  =  108,  M  =  100) 

(e/h  =  27,  M  =  400) 

(e/h  =  54,  M  =  400) 

err_hl  for  (e/h  =  27,  M  =  100) 

9.799  10"2 

9.095  10”2 

8.961  10-2 

(e/h  =  108,  M  =  100) 

(12111111 

(e/h  =  27,  M  =  400) 

(e/h  =  54,  M  =  400) 

err_hl  for  cf’l^1'Vu^’^r  ( e/h  =  27,  M  =  100) 

6.000  10~2 

4.542  10"2 

3.018  10”2 

(e/h  =  108,  M  =  100) 

HKitifelltBa 

(e/h  =  27,  M  =  400) 

'oT 

Oi 

HP- 

o 

o 

Table  1:  Approximation  of  E (Vue)  in  the  L2-norm  (err_hl  )  by 

and  cf/^'Vu^  in  function  of  N  (the  various  lines  correspond  to 
various  values  of  h  and  M ). 
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